{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>*Introduction*</h2>\n",
    "\n",
    "This file enables a user to construct and manipulate geometric objects in $\\mathbb{R}^3$. The constructions and manipulations are performed using a conformal model of $\\mathbb{R}^3$. A user need not know much about the conformal model, as all constructions and manipulations are via the functions provided here.   \n",
    "\n",
    "My intent is that the functions are self documenting through their code and comments. \n",
    "\n",
    "Vectors passed to the functions must be in conformal representation. \n",
    "Exceptions: pt, which converts a 3D point to a conformal point, a 3D normal vector $\\mathbf{n}$, and a 3D parallel bivector $\\mathbf{B}$.\n",
    "\n",
    "Objects returned by the functions are in conformal representation. Exception: tp, which returns a 3D point.\n",
    "\n",
    "To start, pull down the \"Cell\" menu item and choose \"Run All\".\n",
    "\n",
    "Comments and proposed changes or additions are welcome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Conformal Model, Amsterdam convention.  Dorst et al. p. 361\n",
    "from sympy import *\n",
    "from galgebra.ga import Ga\n",
    "from galgebra.mv import *\n",
    "# from lt import *\n",
    "# from sympy import *\n",
    "\n",
    "cm3coords = (o,x,y,z,infty) = symbols('o 1 2 3 infty', real=True)\n",
    "cm3g = '0 0 0 0 -1, 0 1 0 0 0, 0 0 1 0 0, 0 0 0 1 0, -1 0 0 0 0'\n",
    "cm3 = Ga('o e_1 e_2 e_3 oo', g = cm3g, coords = cm3coords)\n",
    "(eo, e1, e2, e3, eoo) = cm3.mv()\n",
    "ep = eo - eoo/2  # ep^2 = +1  GACS 408\n",
    "em = eo + eoo/2  # em^2 = -1\n",
    "E = eo^eoo\n",
    "Ga.dual_mode('Iinv+')\n",
    "#cm3coords = (o,x,y,z,infty) = symbols('o x y z \\infty', real=True)\n",
    "#cm3 = Ga('o e_x e_y e_z \\infty', g = cf3g, coords = cf3coords)\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pt(arg): # R^3 vector --> conformal point. \n",
    "    if isinstance(arg,str):           # Return general 3D point\n",
    "        v = cm3.mv(arg, 'vector')     # General conformal vector \n",
    "        v = v + (v < eoo)*eo + (v < eo)*eoo  # 3D part \n",
    "        v = eo + v + (v<v)*eoo/2\n",
    "    elif arg == 0:\n",
    "        v = eo\n",
    "    elif (arg < eoo) == 0:    # Return point for 3D vector in arg\n",
    "        v = eo + arg + (arg<arg)*eoo/2\n",
    "    else: v = arg     # arg already in conformal representation   \n",
    "    return(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tp(arg): # conformal point --> R^3 vector\n",
    "    if isinstance(arg,str):   # Return general 3D vector\n",
    "        v = cm3.mv(arg, 'vector')\n",
    "    else:                     # Return 3D vector part of arg\n",
    "        v = arg\n",
    "    v = v + (v < eoo)*eo + (v < eo)*eoo\n",
    "    return(v)       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize(v): \n",
    "    if (v < eoo) == 0: # Normalize 3D vector\n",
    "        return(v/sqrt((v<v).scalar()))\n",
    "    else:  # Normalize conformal vector: set eo coeff to 1.\n",
    "        return(-v/(v<eoo))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scalar(arg):\n",
    "    return(cm3.mv(arg, 'scalar')) # Save user from typing all this"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>* Create direct representations of geometric objects *</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def round(*args):  # args are conformal points\n",
    "    ans = args[0]\n",
    "    for i in range(1,len(args)):\n",
    "        ans = ans ^ args[i]\n",
    "    return(ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flat(*args):   # args are conformal points\n",
    "    return(round(*args) ^ eoo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line(p,q): # If q is 3D, line thru p parallel to q returned\n",
    "    return(flat(p,q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plane(p,q,r):\n",
    "    return(flat(p,q,r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def circle(p,q,r):\n",
    "    return(round(p,q,r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sphere(p,q,r,s):\n",
    "    return(round(p,q,r,s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>* Create dual representations of geometric objects *</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dualLine(p, B): # thru point p, orthogonal to 3D bivector B\n",
    "    return(p < (B*eoo))  # A vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dualPlane(p,n):    # n: GA^3 normal vector    \n",
    "    m = normalize(n)\n",
    "    if isinstance(p,(int, long, float)):\n",
    "        p = scalar(p)         # Python scalar -> GAlgebra scalar\n",
    "    if (p!=0) and ((p<p)==0): # p: point on plane. \n",
    "        return(p < (m^eoo))   # a vector\n",
    "    else:                     # p: distance to origin.\n",
    "        return(m + (p*eoo))   # a vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dualSphere(c,rho):  # c:center. \n",
    "    if isinstance(rho,(int, long, float)):\n",
    "        rho = scalar(rho)   # Python scalar -> GAlgebra scalar\n",
    "    if (rho!=0) and ((rho<rho)==0):  # rho: point on sphere \n",
    "        return(rho < (c ^ eoo))  \n",
    "    else:                            # rho: radius. \n",
    "        return(c - (rho*rho*eoo)/2)  # A vector     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dualCircle(c,rho,n): # c:center. rho:radius. n:normal vector\n",
    "    ds = dualSphere(c,rho)\n",
    "    dp = dualPlane(c,n)    \n",
    "    return(ds^dp)          # A BIvector "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>*  Geometric operations *</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def translate(object,a3): # a3: 3D vector\n",
    "    return(1 - a3*eoo/2)*object*(1 + a3*eoo/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotate(object,itheta):\n",
    "    return(exp(-itheta/2)*object*exp(itheta/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def invert(p, norm=False):   # GACS 513\n",
    "    ans = -(eo - eoo/2)*p*(eo - eoo/2) \n",
    "    if norm:\n",
    "        ans = normalize(ans)\n",
    "    return(ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reflect point p in hyperplane with normal 3D vector n.\n",
    "def reflect(p,n):\n",
    "    return(-n*p*(n/norm2(n)))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Can be considerably simplified: A Covariant Approach ..., 16 \n",
    "def dilate(p, alpha, norm = False):  # Dilate by alpha (> 0)\n",
    "    ans = exp(E*ln(alpha)/2)*p*exp(-E*ln(alpha)/2)\n",
    "    if norm:\n",
    "        ans = normalize(ans)\n",
    "    return(ans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>* Play *</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
